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Abstract 

Background: Asthma is a chronic inflammatory airway disease influenced by genetic and environmental factors that 
affects -300 million people worldwide, leading to -250,000 deaths annually. Glucocorticoids (GCs) are well-known 
therapeutics that are used extensively to suppress airway inflammation in asthmatics. The airway epithelium plays an 
important role in the initiation and modulation of the inflammatory response. While the role of GCs in disease 
management is well understood, few studies have examined the holistic effects on the airway epithelium. 

Methods: Gene expression data were used to generate a co-transcriptional network, which was interrogated to 
identify modules of functionally related genes. In parallel, expression data were mapped to the human protein-protein 
interaction (PPI) network in order to identify modules with differentially expressed genes. A common pathways 
approach was applied to highlight genes and pathways functionally relevant and significantly altered following GC 
treatment. 

Results: Co-transcriptional network analysis identified pathways involved in inflammatory processes in the epithelium 
of asthmatics, including the Toll-like receptor (TLR) and PPAR signaling pathways. Analysis of the PPI network identified 
RXRA, PPARGC1A, STAT1 and IRF9, among others genes, as differentially expressed. Common pathways analysis 
highlighted TLR and PPAR signaling pathways, providing a link between general inflammatory processes and the 
actions of GCs. Promoter analysis identified genes regulated by the glucocorticoid receptor (GCR) and PPAR pathways 
as well as highlighted the interferon pathway as a target of GCs. 

Conclusions: Network analyses identified known genes and pathways associated with inflammatory processes in the 
airway epithelium of asthmatics. This workflow illustrated a hypothesis generating experimental design that integrated 
multiple analysis methods to produce a weight-of-evidence based approach upon which future focused studies can be 
designed. In this case, results suggested a mechanism whereby GCs repress TLR-mediated interferon production via 
upregulation of the PPAR signaling pathway. These results highlight the role of interferons in asthma and their 
potential as targets of future therapeutic efforts. 
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Background 

Asthma is a chronic inflammatory airway disease influ- 
enced by genetic and environmental factors, character- 
ized by variable and recurring symptoms (wheezing, 
coughing, chest tightness, shortness of breath), airflow 
obstruction and bronchospasm [1]. Asthma affects -300 
million people worldwide, leading to -250,000 deaths 
annually [2]. The majority of the asthmatic population is 
characterized by the development of an inflammatory re- 
sponse following allergen exposure (allergic asthma). 
The airway epithelium plays a critical role in allergen 
sensitization in asthma. Allergens and allergen products 
activate specific receptors in the surface of the epithelial 
cells, including Toll-like receptors (TLRs). The activated 
epithelium releases inflammatory cytokines that enhance 
dendritic cell migration to lymph nodes, recruitment of 
innate immune and Th2 cells to the airways, and activa- 
tion of Th2 cells [3,4]. 

While no cure exists to date, glucocorticoids (GCs) are 
commonly prescribed to control asthma symptoms [5,6]. 
Glucocorticoids bind to the glucocorticoid receptor 
(GCR), which interacts and represses (trans-repression) 
the activity of transcription factors that regulate the ex- 
pression of inflammatory genes {e.g., transcription factors 
API and NFKB1). The activated GCR also binds to gluco- 
corticoid response elements (GRE) in the DNA and pro- 
motes the expression of genes with anti-inflammatory 
actions {e.g., GILZ, MKP-1) [7]. 

Although the effects of GCs on the airway epithelium 
are well known, the mechanisms mediating GCR actions 
are not as well understood. A significant percentage of 
the asthmatic population requires the maximum doses 
of inhaled GCs and suffers the associated side effects 
and a small fraction of asthmatics are resistant to GCs 
[8]. Accordingly, increased insight into the molecular 
mechanisms mediating GC anti- inflammatory actions 
will assist in the development of therapies that minimize 
resistance and side effects as well as increase our under- 
standing of the disease. In addition, an understanding of 
individual responses to GC treatment will assist in estab- 
lishing treatment regimens tailored for personal suscep- 
tibility to GCs. 

We previously investigated the effects of a synthetic 
GC, fluticasone propionate (Flovent), on the transcrip- 
tional profile of the airway epithelium of asthmatics [9]. 
In that study the aim was to identify genes regulated by 
GCs in the airways, and therefore the focus was on the 
most significantly altered genes. We demonstrated the 
GC-dependent regulation of a number of genes, includ- 
ing CLCA1, POSTN, SERPINB2 and FKBPS [9]. To gain 
further insight into the processes governing the actions 
of GCs in the airways, we performed an unbiased ana- 
lysis to interrogate the microarray dataset from a 
systems-based perspective. Our approach was primarily 



that of a hypothesis generating experiment, where two 
complementary network approaches were applied to mine 
different aspects of cellular networks. The first method 
used co-transcriptional information to characterize the 
functional pathways activated in the airways. The second 
method combined protein-protein interaction (PPI) net- 
works with gene expression changes to identify pathways 
altered following Flovent treatment. Subsequently, com- 
bining results from both networks provided a weight of 
evidence approach to elucidate mechanistic processes in 
the action of Flovent in asthmatics. These results were fur- 
ther supported by promoter enrichment analysis to 
highlight a putative role of PPAR in suppressing the TLR- 
mediated interferon pathway following GC treatment. The 
combination of these different approaches demonstrated 
how multiple analysis methods can be combined in an in 
silico experiment to generate novel hypotheses. Taken to- 
gether, these efforts provided a concerted investigation 
into the action of GCs in the airway epithelium of 
asthmatics. 

Results 

Statistical analysis 

Figure 1 shows an overview of the analysis workflow fol- 
lowed in this study. Linear model analysis of differences 
due to Flovent treatment revealed 1,492 genes that were 
differentially expressed in Flovent versus placebo treated 
subjects {P<0.05, Additional file 1: Table SI). The P- 
values were not corrected for multiple testing and there- 
fore the list of genes is referred to as 'RAW'. This list 
included many genes with a purported role in regulation 
of the immune response in asthma: interleukins {IL2, 
IL27 and IL33), chemokines {CXCL1, CXCL17, CCLS 
and CCL27), mucins {MUC2, MUC13 and MUCL1), 
TLR7 and NOS2, among others. Functional enrichment 
analysis resulted in 13 significantly enriched KEGG 
pathways and 132 GO biological processes {P<0.05) 
(Additional file 2: Table S2). Additional file 3: Figure SI 
shows the 10 KEGG pathways that resulted after filtering 
(see Methods), including the top ones Allograft rejec- 
tion', Antigen processing and presentation' and 'Glutathi- 
one metabolism'. Results for GO enrichment included: 
antigen processing and presentation of peptide antigen 
via MHC class I', antigen processing and presentation' 
and positive regulation of B cell proliferation' as the top 
ranking biological processes. Although the list of uncor- 
rected P-values (RAW) contained many genes and path- 
ways related to asthma, the results possessed a potential 
high false positive rate (16,015 x 0.05 = 800 theoretical 
false positives). After correction for multiple testing by 
the methods of Benjamini and Hochberg [10], 115 genes 
remained significantly altered with q<0.05 (Additional 
file 1: Table SI). This list of genes is referred to as 'FDR', 
and subsequent KEGG pathway enrichment resulted in 
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Figure 1 Overview of the analysis workflow. This figure shows an overview of the analysis workflow followed in this study as described in the 
text. From the initial microarray data (orange), three different approaches were pursued (light green). (I) Statistical analysis (II) Co-transcriptional 
network analysis and (III) Interactome network analysis. Each of these paths ends with a functional analysis (blue), which was combined into a 
common pathways analysis at the end of the pipeline (purple). Promoter analysis (lavender) was performed afterwards as a validation step for the 
conclusions from the common pathways analysis. LM = linear models. RAW = unadjusted P-values (not corrected for multiple hypothesis testing). 



five pathways overrepresented; 'Steroid hormone biosyn- 
thesis', 'Sulfur metabolism', 'Primary bile acid biosynthesis', 
'Endocytosis' and 'Retinol metabolism' (Additional file 2: 
Table S2 and Additional file 3: Figure SI). GO enrichment 
analysis returned 123 significant terms including negative 
regulation of B cell activation' and 'B cell mediated im- 
munity' and other biological processes related to asthma. 
However, only a few genes were associated with these 
selected terms (1-2 genes with P~ 1E-03). 

Analysis of the co-expression network 

A co-transcriptional network was constructed from the 
microarray data using the ARACNE method [11]. The net- 
work was built based upon within subject differences in 
transcript levels between the first sample (pre-treatment) 
and the second sample (after Flovent treatment or placebo). 
The inferred network contained 7,676 nodes and 49,499 
edges (Additional file 3: Figure S2), and the degree distribu- 
tion followed a power-law with y = 1.69 (Additional file 3: 
Figure S3). The MCODE method was applied to identify 



clusters of genes with similar expression profiles in the in- 
ferred network [12]. This method was originally designed 
to detect complexes in PPI networks, but can also be used 
to detect clusters in co-transcriptional networks. MCODE 
retrieved 109 modules, with module 1 (Ml) being the 
highest ranked {i.e., greatest connectivity in the network). 
Module Ml contained 150 highly correlated genes 
(average Spearman p = 0.8074; Figure 2). Pathway analysis 
resulted in 19 significantly enriched KEGG pathways 
(P< 0.05, Additional file 2: Table S2), which were reduced 
to 11 after pathway filtering (Additional file 3: Figure SI). 
Most of the pathways enriched in Ml are directly related 
to biological processes in asthma, including Asthma', 'Cell 
adhesion molecules (CAMs)! 'Antigen processing and 
presentation! 'Cytokine-cytokine receptor interaction! 
'Toll-like receptor signaling pathway! 'PPAR signaling path- 
way! 'Leukocyte transendothelial migration' and 'Chemokine 
signaling pathway'. GO enrichment analysis revealed 313 
significantly enriched biological processes (P<0.05), with 
many categories related to asthma pathogenesis located 
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GeneOntology terms 

□ activation of immune response 

□ activation of plasma proteins involved in acute 
inflammatory response 

□ complement activation, classical pathway 

□ humoral immune response 

□ immune response 

□ inflammatory response 

■ protein processing 

■ regulation of response to stimulus 

□ response to external stimulus 

■ response to stress 

KEGG Pathways 

□ Antigen processing and presentation 

□ Asthma 

□ Cell adhesion molecules (CAMs) 

□ Chemokine signaling pathway 

□ Cytokine-cytokine receptor interaction 

□ Intestinal immune network for IgA production 

■ Leukocyte transendothelial migration 

□ Lysosome 

■ PPAR signaling pathway 

■ Toll-like receptor signaling pathway 

■ Type I diabetes mellitus 

Figure 2 MCODE module Ml. Selected KEGG and GO enriched terms in the MCODE module Ml. This module contains the greatest number of 
genes of all modules identified via MCODE from the co-transcriptional network and has the highest MCODE score {i.e., connectivity density). A 
red label indicates a gene annotated in GO (biological process) with transcription factor activity. This module contains no differentially expressed 
genes identified at P<0.05. 



top in the ranking (i.e., had lowest P-values). In particular, 
the most significant term was 'inflammatory response' 
with 28 associated genes (P= 1.12xl0~ 19 ), including che- 
mokine ligands (CXCL3, CCL4, CCL18), chemokine re- 
ceptor (CCR1), complement components (C7QA, C1QB, 
C1QC, C2), integrins (ITGB2), interleukins (IL1B) and 
lymphocyte antigens (LY86). Other top-ranked terms 
included 'leukocyte migration; regulation of T cell activa- 
tion', neutrophil chemotaxis' and regulation of interleukin 
2 production (Additional file 2: Table S2). 

The remaining modules generally contained less infor- 
mation, with the exception of modules M5, M6 and 
M9. Module M5 contained the greatest number of 
differentially expressed genes (Additional file 3: Figure 
S4, 26 differentially expressed genes at P<0.05, and 22 
genes with q < 0.05), which were both negatively and 
positively correlated (absolute mean Spearman p = 0.665; 
Additional file 3: Figure S5). Only three KEGG pathways 
were overrepresented, including 'Steroid hormone biosyn- 
thesis' (Additional file 2: Table S2). GO functional analysis 
revealed 40 biological processes; however, 35 of them con- 
sisted of only a single gene association. Overrepresented 
terms included steroid metabolic process', 'thiamin trans- 
port; 'positive regulation of epithelial cell proliferation' and 
'progesterone receptor signaling pathway' (Additional file 2: 
Table S2). Module M6 contained 39 correlated genes (aver- 
age Spearman p = 0.6730; Additional file 3: Figures S5 and 
S6) and appeared to be functionally related to module Ml. 
Some of the genes in module M6 participate in pathways 
enriched in module Ml, including TLR4 ('Toll-like receptor 
signaling pathway), PPARG and NFAM1 ('PPAR signaling 



pathway'). Module M9 had 7 correlated genes (average 
Spearman p = 0.8055) genes, including JUN, JUNB and 
FOS, which are part of the API transcription factor that 
regulates the expression of inflammatory genes, and is a 
target of the activated GCR. The enriched KEGG path- 
ways that passed the filtering were 'B-cell receptor signaling 
pathway! 'MAPK signaling pathway! 'T cell receptor signal- 
ing pathway' and 'Toll-like receptor signaling pathway' 
(Additional file 3: Figure S7; Additional file 2: Table S2). 
The top-ten ranked GO terms included 'response to cyto- 
kine stimulus' and 'SMAD protein signal transduction' 
(Additional file 2: Table S2). 

Analysis of the human interactome 

Protein-protein interaction (PPI) networks represent the 
skeleton on which signaling pathways are constructed and 
can be used to identify modules of differentially expressed 
genes related to the same or overlapping signaling path- 
ways [13]. The human interactome was obtained from the 
BioGRID database, and the BioNet software was used to 
identify modules significantly enriched in differentially 
expressed genes [14,15]. BioNet detected a subnetwork 
containing 32 genes, including 19 upregulated and 4 
downregulated {P<0.05, 19 with q<0.05) by Flovent 
treatment (Figure 3). KEGG pathway enrichment analysis 
of this module highlighted 9 significant pathways includ- 
ing 'Toll-like receptor signaling pathway! 'Steroid hormone 
biosynthesis; 'PPAR signaling pathway! 'ErbB signaling path- 
way! and 'Adipocytokine signaling pathway' (Additional file 
3: Figure SI; Additional file 2: Table S2). GO analysis identi- 
fied 207 processes, with 'alpha-beta T cell activation during 
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GeneOntology terms 

n cellular response to hormone stimulus 

□ negative regulation of mast cell cytokine production 

□ negative regulation of T-helper 2 cell differentiation 

□ negative regulation of T-helper 2 type immune response 

□ regulation of cytokine production during immune response 

□ regulation of inflammatory response 

■ regulation of T-helper cell differentiation 

■ response to cytokine stimulus 

■ response to peptide hormone stimulus 

KEGG Pathways 
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□ ErbB signaling pathway 

□ PPAR signaling pathway 

■ Steroid hormone biosynthesis 
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Figure 3 BioNet module. BioNet module annotated with selected KEGG pathways and GO terms obtained from functional enrichment analysis. 
Differentially expressed genes (P<0.05) are indicated with a red circle (up-regulated in response to Flovent treatment) or a blue circle (down- 
regulated in response to Flovent treatment). A red label indicates a gene annotated in GO (biological process) with transcription factor activity. 
Grey circles indicate genes not expressed (based on XIST filtering - see Methods), which are connected to other genes through gray dashed lines. 



immune response' as the top-third term, although it only 
had a single associated gene. Some other significant bio- 
logical processes not in the top rank have a direct relation 
to asthma pathogenesis, including negative regulation of T- 
helper 2 type immune response^ response to cytokine 
stimulus,' regulation of inflammatory response,' 'SMAD pro- 
tein signaling transduction' and cellular response to hor- 
mone stimulus' (Additional file 2: Table S2). However, these 
terms also show a low number (1-2) of associated genes. 

Analysis of common pathways 

Networks were constructed using the two different 
approaches of co-transcriptional and PPI-based analyses. 
To identify cellular processes active in the airways of 
asthmatics that are targeted by GC treatment, we exam- 
ined common pathways found in the different modules. 
We focused on common pathways between the highest 
ranked MCODE module (Ml) and the PPI module 
(BIONET). The pathways, along with the number of 
genes in each pathway, are shown in the heatmap in 
Figure 4. The only common pathways were 'Toll-like re- 
ceptor signaling pathway' and 'PPAR signaling pathway'. 
To determine whether these two pathways appeared in 
common by chance, we performed a simulation with 
random lists of genes of the same size as the BIONET 
and Ml modules. For each random list, we performed 
enrichment analysis and determined the common path- 
ways. By repeating this process 1,000 times and comput- 
ing the number of common pathways similar to those 
detected in our analysis of the Flovent dataset, we esti- 
mated the probability of observing common pathways by 
chance. The simulation resulted in a probability of ran- 
dom appearance as common pathways of P = 0.002 for 
'Toll-like receptor signaling pathway'. 'PPAR signaling 



pathway' did not appear as a common pathway in any of 
the simulations {P=0). These results demonstrate the 
robustness of the common pathways approach. 

Promoter enrichment analysis 

The analysis of common pathways highlighted the TLR and 
PPAR pathways as targets of GC actions. Consequently, 
some of the genes regulated by these pathways should be 
differentially expressed following Flovent treatment. To ex- 
plore this possibility we analyzed the promoter region of 
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Figure 4 Common pathways analysis. Analysis of common 
pathways highlighted pathways present in the BioNet and selected 
MCODE modules. This approach resulted in Toll-like receptor 
signaling and PPAR signaling pathways as the pathways in common 
with other modules. See Additional file 3: Figure S1 for a full 
summary of the functional enrichment analysis. 



Diez et al. BMC Medical Genomics 201 2, 5:27 
http://www.biomedcentral.eom/1755-8794/5/27 



Page 6 of 13 



differentially expressed genes (q<0.05) to look for tran- 
scription factor binding sites (TFBS) associated with TLR 
and PPAR pathways. The method Pscan was used to find 
enriched TFBS from the JASPAR and TRANSFAC data- 
bases [16]. The top enriched motifs in the JASPAR database 
were IRF2 and IRF1 (P = 4.1E-4 and P = 4.1E-3), whereas 
the top motif in the TRANSFAC database was V$IRSE_01 
(P = 4.1E-3) (Additional file 4: Table S4). Motifs V$IRF1 
(P = 0.027), V$IRF2 (P = 0.041) and V$IRF7 (P = 0.07) were 
also significantly enriched in the TRANSFAC database. 
ISRE (Interferon Stimulated Response Element) motifs are 
involved in the regulation of gene expression in response to 
interferon type 1 (IFNA-alpha/beta) activation [17]. IRF 
motifs are the target of Interferon Regulatory Factors 
(IRFs), which function as transcriptional co-activators of 
interferon-regulated genes [18,19]. Interestingly, type 1 
interferons (alpha and beta) are produced after activation of 
TLRs (especially TLR4 and TLR7). The TFBS results also 
showed enrichment of PPARG/RXRA binding motifs from 
the JASPAR database (PPARG:::RXRA, P = 0.0182) and the 
TRANSFAC database (V$PPARG_01; P = 0.084). 

The enrichment motifs determined by Pscan provided 
the possibility that a transcription factor could bind to the 
promoter of genes. Cscan (http://159.149.109.9/cscan/) 
was used in order to obtain some experimental evidence 
for binding of transcription factors to the promoter region 
of the list of differentially expressed genes [20]. Cscan 
compiles ChlP-seq data from published experiments 
assessing the binding of several transcription factors and 
epigenetic markers in different cell lines for human and 
mice. Results from Cscan identified GCR as the most sig- 
nificantly enriched transcription factor in six different cell 
lines subjected to variable doses of a glucocorticoid recep- 
tor agonist (Additional file 4: Table S4). STAT1 and RXRA 
binding were detected at the respective ranks of 10 and 12 
(Bonferroni corrected P- values: q = 0.0009 and q = 0.0018, 
respectively). Binding for PPARs or IFR9 could not be 
assessed due to a lack of available experimental data. 

Interferon regulatory network 

Promoter analysis suggested the type 1 interferon pathway 
as the most important target of GC actions. To further 
validate this finding we queried the database Interferome 
(http://www.interferome.org), which maintains a curated 
list of known interferon-regulated genes extracted from 
studies on which control and interferon-treated samples 
were analyzed. We found that 20 out of 114 genes in the 
FDR list (17%) have been previously identified as being 
regulated by interferons (Additional file 5: Table S5). This 
results in an associated probability of rejecting the null 
(no enrichment) hypothesis based on the hypergeometric 
distribution of P = 0.0505, indicating that the FDR list is 
enriched in interferon-regulated genes. Additional file 5: 
Table S5 contains the number of datasets in which 



regulation by a certain type of interferon was found. Type 
1 interferons regulated 17 out of 20 genes in a total of 72 
datasets, whereas type 2 interferons were the only source 
of regulation for three genes in a total of 3 datasets. Taken 
together, these results suggest that GCs target the type 1 
interferon pathway in the epithelium of asthmatics. 

Discussion 

Classical statistical and functional enrichment analysis of 
the Flovent/placebo dataset provided little insight into 
the biological processes associated with the treated asth- 
matics. Results from the uncorrected tests contained 
many genes known to be associated with asthma, and 
some enriched pathways and GO terms linked these 
genes to inflammatory responses and metabolism of 
GCs. However, the high number of potential false posi- 
tives hampered the results and correction for multiple 
testing eliminated most of these genes from the final list. 
Genes selected after multiple testing correction (FDR) 
identified the SERPINB2 and FKBP5 genes, which were 
previously validated [9]. However, the FDR approach 
failed to identify CLCA1, another gene experimentally 
confirmed [9], although it was significantly altered based 
on uncorrected P-values. This effect could partially be 
attributed to the reduced statistical power associated 
with the use of multiple testing correction methods 
compared to the uncorrected analysis [21]. It has been 
well demonstrated that the excessive penalty imposed on 
P-values following FDR adjustment increases the num- 
ber of false negatives [22]. On the other hand, POSTN, 
which was also experimentally validated, had a P>0.05 
and therefore was not differentially expressed in our 
analysis. This comparison highlights some of the limita- 
tions associated with statistical analyses that rely on a 
single source of evidence to derive conclusions on bio- 
logical function, and argues for the need for integrative 
approaches to data analysis. 

Network analysis offers some advantages over classical 
analysis by being able to incorporate additional informa- 
tion from multiple sources [21,23,24]. We applied a 
weight-of-evidence approach by using multiple methods 
that mine different aspects of the cellular processes, 
which converged to 'Toll-like receptor signaling path- 
way and 'PPAR signaling pathway as the most relevant 
pathways related to GC treatment (Figure 4). Toll-like 
receptors (TLRs) regulate inflammatory responses by 
inducing the expression of inflammatory cytokines upon 
binding of viral or bacterial proteins and mediate the 
signaling pathways that regulate innate and Th2 
responses in the epithelium [3,25,26]. Toll-like receptors, 
including TLR4 (found in module M6) and TLR7 (found 
in module Ml) mediate the production of interferons 
alpha and beta (type 1 interferons) [27]. This important 
link between TLR and the interferon 1 pathway, not 
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evident from the pathway enrichment analysis, was revealed 
by the promoter enrichment analysis. 

Type 1 interferons act in an autocrine or paracrine 
manner and bind to interferon-alpha receptors in the 
membrane. This activates the Janus kinase and signal 
transducers and activators of transcription (JAK/STAT) 
pathway, which regulates the expression of inflammatory 
genes [17]. Phosphorylation of STAT1 by JAK1 activates 
STAT1, which forms a heterodimer with STAT2 and 
recruits IRF9 to form the IFN-stimulated gene factor 3 
(ISGF3) transcriptional complex [18,19]. The analysis of 
the human interactome identified STAT1 and IRF9, 
which are downregulated (P< 0.05) by Flovent (Figure 3). 
In addition, promoter analysis identified ISRE as the top 
enriched transcription factor motif in the list of differen- 
tially expressed genes (q<0.05). This suggests a mech- 
anism in which GCs regulate the type 1 interferon 
pathway by downregulating the proteins involved in 
interferon signal transduction. However, although both 
STAT1 and IRF9 are suggested as differentially expressed 
(P<0.05) and exhibited correlated expression profiles 
(r = 0.771, P = 3.7E-07 > ), their correlation with treatment 
is less clear (Figure 5). Interestingly, GCs are known to 
regulate the type 1 interferon pathway by suppressing 
the phosphorylation of STAT1 [28], providing an alterna- 
tive mechanism of regulation. 

STAT1 phosphorylation also mediates the transduction 
of type 2 interferon (IFN-gamma) signaling by binding 
to gamma-activated sequence (GAS) elements as an 
homodimer. Therefore, it is possible that downregula- 
tion of STAT1 and inhibition of phosphorylation of the 
protein represses type 2 interferon signaling. Indeed, 
GCs are known to inhibit IFN-gamma signaling by 
downregulating STAT1 mRNA and protein expression in 
PBMCs [29]. Gene expression for type 2 interferon and 
its receptor (IFN-gamma receptor, IFNGR) were found 
in the microarray dataset after filtering non-expressed 
genes (not shown). However, promoter analysis only 
identified enrichment of ISRE elements, the binding 



motif for the ISGF3 transcriptional complex, GAS ele- 
ments were not detected. In addition, interrogation of 
the interferome database resulted in a list of interferon- 
regulated genes that were, in the majority of datasets 
assayed, regulated by type 1 interferons. Consequently, 
although a role of type 2 interferons cannot be excluded, 
our findings suggest that type 1 interferons are the target 
of GCs in the epithelium of asthmatics. 

Peroxisome proliferator- activated receptors (PPAR) are 
lipid-activated transcription factors that regulate the ex- 
pression of target genes [30]. PPAR-alpha and -gamma 
modulate allergic inflammation, and agonists are able to 
reduce levels of inflammatory cytokines [31,32]. PPARs 
exert their activity by forming a heterodimer with the 
retinoid receptor RXRA, which then binds to co-activator 
proteins, including PPARGC1A (PGC-i), to regulate gene 
expression. Transcriptional co-activators amplify the tran- 
scription of nuclear receptor regulated target genes [33]. 
In particular, PPARGC1A can recruit other co-activator 
proteins with histone acetyltransferase (HAT) activity that 
open up the chromatin and enhance the expression of tar- 
get genes [34]. Consequently, an increase in the expres- 
sion of RXRA and PPARGC1A co-activators can intensify 
the activity of the PPAR pathway, resulting in an increased 
repression of STAT1 phosphorylation. Both RXRA and 
PPARGC1A are present in the BioNet module and are 
upregulated by Flovent (Figure 3). Promoter analysis also 
resulted in an enrichment of genes with PPARG:RXRA 
motifs (P -0.01). This suggests that GCs modulates the 
activity of the PPAR pathway by upregulating co- 
transcriptional activators. Supporting this evidence, the 
expression profiles for RXRA and PPARGC1A are corre- 
lated (r = 0.4078, P = 0.022) and show a trend of higher 
levels in Flovent-treated patients (Additional file 3: Figure 
S8). 

Our results link GCs to type 1 interferon and PPAR 
pathways. There is previous evidence that connects PPARs 
to interferons, providing an interesting mechanism integrat- 
ing GC actions. Activation of the PPAR-alpha (PPARA) 
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Figure 5 STAT1 and IRF9 expression profiles. STAT1 and IRF9 show correlated expression profiles and are downregulated following Flovent 
treatment. 
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pathway was found to suppress STAT1 phosphorylation in 
rat glia [35]. However, in our study promoter analysis iden- 
tified PPARG, but not PPARA (P ~ 0.2), motifs enriched in 
the list of differentially expressed genes. This lack of signifi- 
cance could be explained by the fact that PPARA motifs 
are missing in the JASPAR database, which retrieved higher 
significant results for the complex PPARG:RXRA transcrip- 
tion factor motif. In addition, our search strategy looked for 
motifs in the 1 kb region upstream of the transcription start 
site, whereas some transcription factors can bind to more 
distal locations. Alternatively, as our data suggest, GC ac- 
tion on type 1 interferons may be mediated via a PPAR- 
gamma-dependent process. For example, it has been shown 
that PPAR-gamma can repress the type 1 interferon path- 
way by downregulating the production of INF-beta upon 
TLR4 activation [36]. Treatment with the PPAR-gamma 
agonist troglitazone and challenge with LPS and poly (I: C) 
impaired IRF3 binding to the IFN-beta-promoter. Downre- 
gulation of IFN-beta prevented activation of the IFN-beta 
receptor and subsequent STAT1 phosphorylation and ISRE 
activation [36]. However, in our dataset type 1 interferons 
were present, but not differentially expressed (average 
P-0.58 for type 1 interferon genes in the microarray). 
Interestingly, activation of the PPAR-gamma pathway was 
previously found to downregulate the expression of IFN- 
gamma activated genes [37]. These findings highlight the 
complex nature of PPAR-mediated interferon regulation, 
which can affect different pathways (type 1 vs. type 2) at 
multiple regulatory points. 

Our findings provide a link between interferon, PPARs 
and GCs that suggests a model similar to that presented in 
Figure 6. In this model, allergens activate the TLR signaling, 
which in turn activates the production of type 1 interferons. 
The alpha/beta interferons then bind to interferon recep- 
tors (IFNAR1/2) that stimulate the phosphorylation of 
STAT1 and promote the expression of genes with a num- 
ber of inflammatory effects. Treatment with GCs upregu- 
lates PPARGC1A and RXRA coactivator molecules, which 
consequently enhance the PPAR pathway. Activation of the 
PPAR pathway inhibits phosphorylation of STAT1 and 
therefore inhibits the interferon pathway. In addition, GCs 
may repress the interferon pathway by downregulation of 
STAT1 and IRF9 transcription factors. In this model, PPARs 
could be potential mediators of the anti-inflammatory 
actions of GCs. Both PPAR-alpha and -gamma inhibit air- 
way inflammation in a murine module of asthma [32], The 
use of PPAR-gamma agonists has been shown to evidence 
improvements in lung function of smokers with asthma 
(improved FEVx; forced expiratory volume in 1 sec), who 
had previously demonstrated a reduced response to GC 
treatment [38]. While the mechanism for the observed 
improvements in lung function was unclear, it was postu- 
lated that PPAR-gamma could independently modulate a 
set of inflammatory genes relative to GCs [39]. 



A link between GCs and the interferon pathways 
has been previously reported [5,40-42], and interfer- 
ons are known to affect symptoms in asthmatics [43- 
46]. Interferon-alpha has been associated with severe 
exacerbation of asthma symptoms [43], which provides a 
simple mechanistic interpretation of the beneficial effects 
of GC-mediated repression of the interferon pathway. On 
the other hand, low doses of interferon-alpha have been 
associated with therapeutic effects in GC-resistant 
patients [44-46]. It is unclear if these disparate responses 
are due to differences in interferon dose or different pa- 
tient phenotype. This clearly highlights the inherent com- 
plexity of the underlying regulatory networks and the 
need of further studies investigating the mechanisms of 
GCs and their relation with the interferon pathway. 

Conclusions 

Multifactorial diseases such as asthma challenge our ability 
to identify distinct mechanisms involved in disease etiology 
and pathology. The advent of high-throughput approaches 
has increased the amount of information that can be 
extracted from a sample set, but linking together disparate 
datasets remains an elusive goal. Technologies that summa- 
rize the principal disease processes affected {e.g., functional 
enrichment analysis) assist in characterizing the primary 
functional features, but often fail to suggest detailed infor- 
mation of disease and treatment action. Accordingly, ap- 
proaches that incorporate information regarding multiple 
biological aspects from disparate sources can stimulate the 
generation of hypotheses that will uncover novel mechan- 
isms. Given that asthma is a heterogeneous inflammatory 
condition, it is of interest to examine for distinct molecular 
phenotypes. For example, it has been shown that respon- 
siveness to inhaled GCs correlates with the degree of Th2 
inflammation. It is possible that distinct subgroups of asth- 
matics exist with specific shifts in PPAR/TLR pathways that 
correlate with GC responsiveness or other molecular phe- 
notypes, such as the IL-13 related Th2 sub-phenotypes pre- 
viously described for the current data set [47,48]. Further 
investigations to investigate sub-phenotypical responses in 
the PPAR/TLR pathways are thus warranted. 

With the workflow presented herein, we were able to 
identify specific pathways as being targets of GCs in the 
epithelium of asthmatics. Results highlighted the promin- 
ent role of interferons in mediating inflammatory processes 
in asthmatics, and further supported the finding that GCs 
mediate their activities by repressing the interferon path- 
way. We found that this repression may be mediated by 
GC-dependent upregulation of the PPAR pathway. A po- 
tential mechanism of direct repression of the interferon 
pathway by downregulation of key transcription factors is 
also suggested. If confirmed, these findings will be poten- 
tially valuable in the design of new therapies tailored for 
the GC-resistant asthmatic sub-population. 
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Figure 6 Model of glucocorticoid (GC) actions in the airways epithelium. Simplified (A) and detailed (B) model of PPAR-mediated anti- 
inflammatory actions of GC in lung epithelial cells. Following activation by allergens, TLRs activate the type 1 interferon pathway. Interferons 
induce the expression of inflammatory genes mediated by phosphorylated STAT1 and IRF9 transcription factors. GC treatment represses the 
interferon pathway via two potential mechanisms. First, upregulation (orange boxes) of RXRA and PPARGC1A enhances the activity of the PPAR 
pathway, which inhibits STAT1 phosphorylation, repressing interferon signaling. Secondly, down-regulation (blue boxes) of STAT 1 and IRF9, reduces 
the amount of the transcriptional complex responsible for interferon-regulated gene expression. Clear arrowheads indicate phosphorylation. Solid 
arrowheads indicate activation. A dashed line represents transcriptional regulation (arrow: activation; T: repression). 



Materials and methods 

Microarray data acquisition and preprocessing 

The microarray data set was published previously, and 
is available in the GEO database (accession number 
GSE4302) [9]. The dataset consisted of 64 Affymetrix 
HGU133plus2 arrays measuring expression profiles from 
airway epithelial brush biopsies obtained from 32 asth- 
matics in a randomized, blinded parallel group study design 
with sampling at baseline levels and after 1 week treatment 
with either placebo or fluticasone propionate (Flovent). Ob- 
servation 55 was identified as an outlier due to poor tech- 
nical quality of the microarray data, and was excluded from 
all analyses. Custom CDFs (Chip Description File) from the 
BrainArray Project were used to map the array probes to 
17,788 probe sets corresponding to individual Entrez Gene 
entries [49]. Arrays were preprocessed with the RMA 
method from the Bioconductor package affy [50,51]. Non- 
expressed genes were eliminated on the basis of XIST 
expression, which is not expressed in males [52]. The distri- 
bution of intensities for XIST in male and females was com- 
puted (Additional file 3: Figure S9), and the background 
threshold estimated as the median value of the distribution 
in males (log2 median intensity = 5.32). Genes were 



eliminated when the intensity in all arrays was below that 
of the determined threshold value. The resulting microarray 
data contained 16,015 remaining genes. Baseline correction 
was performed by computing the log difference of the in- 
tensity values from baseline and post-treated arrays prior to 
statistical analysis. 

Statistical analysis 

Statistical analysis was performed with the Bioconductor 
package limma [53]. A linear model was constructed in- 
cluding drug (i.e., Flovent or placebo), sex and age in the 
design matrix. P values associated with the comparison 
Flovent vs. placebo were computed. Multiple testing cor- 
rection was applied using the method of Benjamini and 
Hochberg [10]. 

Co-transcriptional network 

A co-transcriptional network was constructed from the 
microarray data using the ARACNE method, which uses 
mutual information (MI) to test statistical dependency be- 
tween expression profiles [11]. Due to concerns about the 
effects of low sample numbers in combination with the in- 
herent biological variability upon the statistical dependency, 



Diez et al. BMC Medical Genomics 201 2, 5:27 
http://www.biomedcentral.eom/1755-8794/5/27 



Page 10 of 13 



all individuals in the dataset (Flovent-treated and placebo- 
treated) were included in the co-transcriptional network in- 
ference. To reduce the inter-individual variance from the 
gene expression profiles, expression values were normalized 
to the baseline values for each individual (i.e., samples col- 
lected prior to initialization of treatment regimen). Mutual 
information was computed using a Gaussian Kernel esti- 
mate. Two critical parameters for the ARACNE method 
are the kernel width and the MI threshold, which are nor- 
mally estimated using precomputed calibration curves. 
However, for this project, custom calibration curves were 
constructed using the MATLAB scripts provided in the 
original publication [11]. A calibration curve for kernel 
width was generated (Additional file 3: Figure S10A) and 
the optimal kernel width was determined to be 0.2179931, 
which is 5.3% lower than that estimated by ARACNE 
(0.230167). The calibration curve for the MI threshold 
was computed (Additional file 3: Figure S10B) and used 
during network reconstruction. Next, the optimal P-value 
cutoff for network reconstruction was determined by 
computing networks at different P-values over a suffi- 
ciently wide range. For each network, the number of 
nodes, edges and degree distribution was computed 
(Additional file 3: Figure SUA), and a P-value of 1.0E-8 
was chosen, which resulted in a power law degree distri- 
bution for the network (Additional file 3: Figures S3 and 
SUB). Finally, in order to assess the effect of sample 
ordering upon correlation estimations (MI), a boot- 
strapped network was computed by setting 100 cycles and 
the consensus network was calculated with a Bonferroni 
corrected q< 0.001 [54]. 

Protein-protein interaction network 

Information about protein-protein interactions (PPI) was 
obtained from the BioGRID database [15], which contains 
experimentally validated data from a range of high- 
throughput and low-throughput techniques. The database 
version 3.1.70 contained information on 9,935 genes, 
which formed a network with 50,531 edges. This network 
contained multiple edges joining the same pair of genes 
{i.e., the interaction has several sources of evidence), and 
therefore for module detection the network was simplified 
by using the simplify function in the R package igraph 
(http://igraph.sourceforge.net), to reduce the number of 
edges per pair of genes to one. The resulting simplified 
network contained 33,282 edges. 

Module detection 

For the co-transcriptional network, detection of clusters 
of densely connected nodes was performed using the 
MCODE plugin for the Cytoscape software with default 
parameters [55]. For the BioGRID PPI network, module 
detection was performed using the Bioconductor pack- 
age BioNet [14], which enables the determination of an 



exact solution to the problem of finding connected sub- 
graphs with low P-values. A Binomial Uniform Mixture 
(BUM) model was fitted to the distribution of P- values, 
and scores were derived for each network node at a 
given FDR (q = 0.05). Finally, these scores were used to 
detect modules using the exact Heinz method [56]. 

Functional enrichment analysis 

Functional enrichment analysis for KEGG pathways and 
Gene Ontology (GO) Biological Process (BP) terms 
(downloaded on 2010/09/07 and 2010/09/04 respectively) 
was performed using a hypergeometric test as implemen- 
ted in the Bioconductor package GOstats [57]. For GO, a 
conditional test that considers the dependence structure 
of the GO terms was performed [58]. The lists of probe 
sets from the statistical analysis and modules from the co- 
transcriptional network were converted into lists of Entrez 
Gene identifiers prior to enrichment analysis. Enrichment 
results were considered significant at P < 0.05. 

KEGG pathways were restricted to those involved in bio- 
logical processes. Consequently, disease pathways were dis- 
carded (except the KEGG pathway Asthma). In addition, 
pathways that were supported with less than two genes in 
any module or list of differentially expressed genes were 
discarded. The final list of KEGG pathways was used to fil- 
ter modules in the co-transcriptional network. Modules 
that did not contain any genes annotated in the KEGG 
pathways that passed the filtering were discarded. Pathways 
were grouped using a hierarchical clustering with Pearson 
correlation as a distance measure, and the gene counts dis- 
played as a heatmap. 

Common pathways analysis 

To link the results from the co-transcriptional and PPI 
network analyses, common pathways were selected when 
present simultaneously in the BIONET module and the 
MCODE module Ml from the co-transcriptional network. 
To assess the possibility that the selected common path- 
ways appeared by chance in our dataset, we performed 
simulations with random lists of genes. First, background 
lists of genes were generated by randomly selecting genes 
from the total number of genes (universe) used for enrich- 
ment analysis. The number of genes in each simulated list 
was the same as the target list {i.e., the same number of 
genes as found in the BIONET and MCODE module Ml). 
A KEGG enrichment analysis was then performed for 
each list as described above, and the common pathways 
were computed. Finally, the common pathways in the 
simulation were compared to those obtained with the real 
dataset and the number of positive hits n recorded. This 
process was repeated N= 1000 times and the probability 
of a pathway appearing by chance was computed as the 
number of times it appeared as common in the simulated 
datasets divided by the number of simulations {P = n / N). 
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Promoter enrichment analysis 

Promoter enrichment analysis was performed to iden- 
tify transcription factor binding motifs enriched in the 
promoter regions of sets of genes by using the soft- 
ware Pscan [16]. Entrez Gene identifiers were con- 
verted into REFSEQ ids and these were loaded into 
the Pscan web server. JASPAR and TRANSFAC data- 
bases of transcription factor motifs were selected to 
analyze the region from -950 to +50 bases relative to 
the transcription start site (TSS). Transcription factor 
binding assessment was performed to identify experi- 
mental evidence (based on published ChlP-seq data) of 
binding of transcription factors to the promoter region 
of the differentially expressed genes. Entrez Gene iden- 
tifiers were converted into REFSEQ ids, which were 
loaded into the Cscan web server [20]. A region be- 
tween -1,000 upstream of the TSS and the transcribed 
region of the gene was selected to identify transcrip- 
tion factor binding events. 

Interferon regulated genes 

The database Interferome (http://www.interferome.org) 
was used to identify interferon target genes [59]. This 
database catalogs genes regulated in microarray and 
other experiments where comparisons between control 
and interferon treated samples were performed. Entrez 
Gene identifiers were converted into EnsEMBL gene 
identifiers and loaded into the Interferome web server. 
Enrichment of interferon-regulated genes was computed 
using the hypergeometric test implemented in the R 
software. The parameters were set with 20,900 as total 
number of genes in the human genome (based on data 
from the EnsEMBL database release 64) and 2,000 as the 
number of interferon-regulated genes (as stated in the 
Interferome web site). Since the number of interferon- 
regulated genes is an estimate, the enrichment calcula- 
tion is only approximate. 
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